Pulse detection from head motions in video

ABSTRACT

Heart rates and beat lengths can be extracted from videos by measuring subtle head motion caused by the Newtonian reaction to the influx of blood at each beat. In an embodiment of the present invention, a method tracks features on the head and performs principal component analysis (PCA) to decompose their trajectories into a set of component motions. The method then selects a component that best corresponds to heartbeats based on its temporal frequency spectrum. Finally, the motion projected to this component is analyzed and peaks of the trajectories are identified, which correspond to heartbeats.

RELATED APPLICATION

This application claims the benefit of U.S. Provisional Application No. 61/839,563, filed on Jun. 26, 2013, entitled “Pulse Detection from Head Motions in Video.” The entire teachings of the above application are incorporated herein by reference.

BACKGROUND

Head movements that are related to cardiac activity are small and mixed in with a variety of other involuntary head movements. From a biomechanical standpoint, the head-neck system and the trunk can be considered as a sequence of stacked inverted pendulums. This structure allows the head unconstrained movement in most axes.

SUMMARY

Several sources of involuntary head movement can complicate the isolation of movements attributable to pulsatile activity. One is the pendular oscillatory motion that keeps the head in dynamic equilibrium. (He, D. et al., “A continuous, wearable, and wireless heart monitor using head ballistocardiogram (bcg) and head electrocardiogram (ecg),” Conf Proc IEEE Eng Med Biol Soc. 2011 (2011). (Hereinafter, “He, D., et al.”)). A second source of involuntary head movement is the bobbing caused by respiration.

The heart is a hollow muscle that pumps blood throughout the body by repeated, cyclical contractions. The heart is composed of four chambers: the left ventricle, right ventricle the left atrium, and the right atrium. During the phase of the cardiac cycle known as diastole, the ventricles relax and allow blood to flow into them from the atria. In the next phase known as systole, the ventricles contract and pump blood to the pulmonary artery and aorta. The aorta, in turn, transports blood to the rest of the body. The head and neck receive blood from the aorta via the common carotid arteries, which further divide into the internal and external carotid arteries in the neck. However, from a biomechanical standpoint, the head-neck system and the trunk can be considered as a sequence of stacked inverted pendulums. This structure allows the head unconstrained movement in most axes, making calculations of the system's motion complicated.

It is uncertain how much detected head motion is attributable to the large acceleration of blood in the aorta compared to the localized acceleration of blood in the carotid arteries moving into the head. He, D, measured ballistocardiac head motions with an accelerometer-based device to be on the order of 10 mG

$\left( {0.098\mspace{14mu}\frac{m}{s^{2}}} \right),$ but it is unclear now much of this movement is attributed to aortic and carotid blood flow forces.

During systole, a healthy adult aorta is capable of impelling a volume of 70-80 ml (0.07-0.08 kg) of blood at considerable speed

$\left( {1\mspace{14mu}\frac{m}{s}} \right)$ and acceleration

$\left( {20\mspace{14mu}\frac{m}{s^{2}}} \right).$ An average adult male weighs 70 kg. Using Newton's 2nd and 3rd laws, the approximate acceleration of the body due to the aortic forces (a_(body)) can be derived:

F_(body) = F_(blood) a_(body) * m_(body) = a_(blood) * m_(blood) ${{a_{body}*70\mspace{14mu}{kg}}} = {20\frac{m}{s^{2}}*0.07\mspace{14mu}{kg}}$ ${a_{body}} = {0.02\mspace{14mu}\frac{m}{s^{2}}}$

Thus, the approximate acceleration of the body due to the aortic forces is

$0.02\mspace{14mu}{\frac{m}{s^{2}}.}$

In addition to aortic blood acceleration, there is a smaller force resulting from blood flowing through the carotid arteries into the head. Although there are studies measuring blood flow, velocity and acceleration of blood in the carotid, experimental measurements on the force imparted by carotid blood flow on the head have not been found. According to one study, blood velocity in the common carotid increases from approximately

${0.02\mspace{14mu}\frac{m}{s}\mspace{14mu}{to}\mspace{14mu} 0.11\mspace{14mu}\frac{m}{s}\mspace{14mu}{in}\mspace{14mu} 0.1\mspace{14mu} s},$ or an acceleration of

$0.9\mspace{14mu}{\frac{m}{s^{2}}.}$ The blood mass transported in this period is roughly 13 ml or 0.013 kg. Assuming both carotid arteries are identical and that the head is 5 kg, we use a similar derivation to the one used in Section 2.1.1 to find that the head should accelerate roughly 0.005 m/s₂ assuming independence from the rest of the body.

Pulse rate captures the average number of cardiac cycles over a period of time (e.g., 30 seconds), which is useful primarily for detecting acute problems. There is a growing body of evidence that measuring beat-to-beat variations provides additional in-formation with long-term prognostic value. The most established of these measures is heart rate variability (HRV). HRV measures the variation in the length of individual normal (sinus) heartbeats and provides an indication of the degree to which the sympathetic and parasymathetic nervous systems modulate cardiac activity. To measure HRV, the interarrival times of beats must be accurately measured, which can be determined by locating the “R” peaks in successive beats in an ECG. A lack of sufficient variation when the subject is at rest suggests that the nervous system may not perform well under stress. Patients with decreased HRV are at an increased risk of adverse outcomes such as fatal arrhythmias.

Pulse rate captures the average heart rate over a period of time (e.g., 30 seconds). It is useful primarily for detecting acute problems. A growing body of evidence shows that measuring beat-to-beat variations provides additional information with long-term prognostic value. (Syed, Z. et al., “Computationally generated cardiac biomarkers for risk stratification after acute coronary syndrome,” Sci. Trans. Med., 3(102):102ra95 (2011). (Hereinafter, “Syed, Z. et al.”)) The most established of these measures is heart rate variability (HRV). HRV measures the variation in the length of individual normal (sinus) heartbeats. HRV provides an indication of the degree to which the sympathetic and parasympathetic nervous systems are modulating cardiac activity. To measure HRV, the interarrival times of beats are accurately measured, which can be determined by locating “R” peaks in successive beats in an ECG. A lack of sufficient variation when the subject is at rest suggests that the nervous system may not perform well under stress. Patients with decreased HRV are at an increased risk of adverse outcomes such as fatal arrhythmias.

Heart rate is a critical vital sign for medical diagnosis. Extracting heart rate without physically contacting the patient is desirable, particularly for populations such as premature neonates and the elderly for whom the skin is fragile and damageable by traditional sensors. Furthermore, for the elderly demographic, continuous, or at least frequent, monitoring outside of clinical environments can provide doctors with timely samples, long-term trends, and statistical analyses. Acceptance of such monitoring depends in part on the monitors being non-invasive and non-obtrusive.

In an embodiment of the present invention, subtle head oscillations that accompany the cardiac cycle are exploited to extract information about cardiac activity from videos. In addition to unobtrusively measuring heart rate, the method can extract other clinically useful information about cardiac activity, such as subtle changes in the length of heartbeats that are associated with the health of the autonomic nervous system.

In an embodiment, a method includes selecting a region of a video, tracking features of the selected region of the video, and analyzing the features of the selected region of the video to determine an oscillation rate of a subject shown in the video.

In an embodiment, the oscillation rate is a heart rate or pulse rate of the subject.

In an embodiment, the method includes extracting a component of the tracked features of the subject video. The extracted component can be a vertical component, horizontal component, or both vertical and horizontal components. However, a person of ordinary skill in the art can recognize the extracted component(s) can have a different orientation based on the input video.

In an embodiment, selecting the region of the video includes finding a region of interest including a head of a subject.

In an embodiment, tracking the features includes determining a maximum distance traveled by each point between consecutive frames of the video and discarding features corresponding to points having a distance in consecutive frames of the video exceeding a mode of distribution.

In an embodiment, analyzing the features can include filtering the extracted features temporally. Analyzing the features can further include generating a plurality of components using principal component analysis on the temporally filtered extracted features. Analyzing the features can include selecting one of the components of the plurality of components generated by the principal component analysis. Analyzing the features can further include detecting peaks of the selected component. Analyzing the features can further include determining the oscillation rate based on the detected peaks of the selected component.

In an embodiment, filtering the extracted features temporally can remove frequencies that are outside a predetermined range of pulse rates.

In an embodiment, performing principal component analysis can include decomposing the filtered features into a set of independent source signals.

In an embodiment, selecting the component includes selecting the component having a clearest main frequency.

In an embodiment, a system can include a selection module configured to select a region of a video, a tracking module configured to track features of the selected region of the video, and an analysis module the features of the selected region of the video to determine an oscillation rate of a subject shown in the video.

BRIEF DESCRIPTION OF THE DRAWINGS

The patent or application file contains at least one drawing executed in color. Copies of this patent or patent application publication with color drawing(s) will be provided by the Office upon request and payment of the necessary fee.

The foregoing will be apparent from the following more particular description of example embodiments of the invention, as illustrated in the accompanying drawings in which like reference characters refer to the same parts throughout the different views. The drawings are not necessarily to scale, emphasis instead being placed upon illustrating embodiments of the present invention.

FIG. 1 is a diagram showing the cyclical movement of blood from the heart to the head via the abdominal aorta and the carotid arteries that cause the head to move in a periodic motion.

FIG. 2 is a flow diagram illustrating a high level overview of an embodiment of the present invention.

FIG. 3 displays the two eigenvectors (ϕ₁ and ϕ₂) for two of the subjects.

FIG. 4 is a diagram illustrating histograms of four of the 16 distributions binned at every 0.05 seconds.

FIG. 5 is a diagram illustrating the results of the subjects, plotting σ_(color) based on σ_(motion).

FIG. 6 is a graph illustrating a plot of σ_(motion) against β for all subjects on a logarithmic scale.

FIG. 7 are frames from videos employed by an embodiment of the present invention.

FIG. 8 is a diagram showing a frame from a 30-second video of a newborn recorded in situ at a hospital nursery and corresponding graphs.

FIG. 9 is a block diagram illustrating an example embodiment of the present invention.

DETAILED DESCRIPTION OF THE INVENTION

A description of example embodiments of the invention follows.

In an embodiment of the present invention, a method takes an input video of a person's head and returns a pulse rate as well as a series of beat locations which can be used for the analysis of beat-to-beat variability. Feature tracking extracts the motion of the head. The method then isolates the motion corresponding to the pulse and projects it onto a 1D signal that allows us to extract individual beat boundaries from the peaks of the trajectory. For this, the method uses PCA and selects the component whose temporal power spectrum best matches a pulse. The method projects the trajectories of feature points onto this component and extract the beat locations as local extrema.

FIG. 1 is a diagram 100 showing the cyclical movement of blood from the heart to the head via the abdominal aorta and the carotid arteries 102 that cause the head to move in a periodic motion. Blood can flow from the heart to the head via the carotid arteries on either side of the head. In an embodiment of the present invention, a method detects pulse from this movement by tracking feature points on a person's head, filtering velocities by a temporal frequency band of interest, and using principal component analysis (“PCA”) to find a periodic signal caused by pulse. The method extracts an average pulse rate from this signal by examining its frequency spectrum and obtains precise beat locations with a simple peak detection algorithm.

The method complements the extraction of pulse rate from video via analysis of the subtle color changes in the skin caused by blood circulation (Poh, M. et al., “Non-contact, automated cardiac pulse measurements using video imaging and blind source separation,” Optics Express, 18(10): 10762-10774 (2010)) (Hereinafter, “Poh”); (Verkruysse, W. et al., “Remote plethysmographic imaging using ambient light,” Optics Express, 16(26): 21434-21445 (2008)) (Hereinafter, “Verkruysse”). The methods described in Poh and Verkruysse average pixel values for all channels in the facial region and temporally filter the signals to an appropriate band. The method of Verkruysse uses these signals directly for analysis, while the method of Poh performs ICA to extract a single pulse wave. Both find the frequency of maximal power in the frequency spectrum to provide a pulse estimate.

Philips also produces a commercial application that detects pulse from color changes in real-time (Phillips, Philips vital signs camera, http://www.vitalsignscamera.com, 2011). (Hereinafter, “Phillips”)). These color-based detection schemes require facial skin to be exposed to the camera. In contrast, an embodiment of the present invention is not restricted to a particular view of the head, and is effective even when skin is not visible. Non-invasive pulse estimation using modalities other than video such as thermal imagery (Garbey, M. et al., “Contact-free measurement of cardiac pulse based on the analysis of thermal imagery,” IEEE Trans Biomed Eng, 54(8): 1418-1426 (2007) (Hereinafter, “Garbey”)) and photoplethysmography (measurement of the variations in transmitted or reflected light in the skin. (Wieringa, F. et al., “Contactless multiple wavelength photoplethysmogtaphic imaging: a first step toward SpO₂ camera technology,” Ann. Biomed. Eng., 33(8): 1034-1041 (2005). (Hereinafter, “Wieringa”)).

Body motion in videos can be analyzed in different medical contexts, such as the measurement of respiration rate from chest movement (Tan, K. et al., “Real-time vision based respiration monitoring system. Proceedings of CSNDSP, pages 770-774 (2010). (Hereinafter, “Tan”.); (Phillips), or the monitoring of sleep apnea by recognizing abnormal respiration patterns (Wang, C. et al. “Vision analysis in detecting abnormal breathing activity in application to diagnosis of obstructive sleep apnoea,” Proceedings of IEEE Eng Med Biol Soc, pages 4469-4473 (2006). (Hereinafter, “Wang”)). Motion studies for diseases include identification of gait patterns of patients with Parkinson's disease (Chen, S. et al., “Quantification and recognition of parkinsonian gait from monocular video imaging using kernelbased principal component analysis,” BioMedical Engineering OnLine, #10 (2011). (Hereinafter, “Chen”)), detection of seizures for patients with epilepsy (Pediaditis, M. et al., “Vision-based human motion analysis in epilepsy—methods and challenges,” Proceedings of IEEE ITAB, pages 1-5 (2010). (Hereinafter, “Pediaditis”)) and early prediction of cerebral palsy (Adde, L. et al. “Early prediction of cerebral palsy by computer-based video analysis of general movements: a feasibility study,” Developmental Medicine & Child Neurology, 52: 773-778 (2010). (Hereinafter, “Adde”)). The movements involved in these approaches tend to be larger in amplitude than the involuntary head movements due to the pulse.

Imperceptible motions in video can also be amplified (Wu, H. et al., “Eulerian video magnification for revealing subtle changes in the world,” ACM Trans. Graph. (Proceedings SIGGRAPH 2012), 31(4) (2012). (Hereinafter, “Wu”); Liu, C. et al., “Motion magnification,” ACM Trans. Graph., 24(3): 519-526 (2005) (Hereinafter, “Liu”)). See also Wu et al., “Linear-Based Eulerian Motion Modulation,” U.S. patent application Ser. No. 13/850,717, Wadhwa et al., “Complex-Valued Phase-Based Eulerian Motion Modulation,” U.S. application Ser. No. 13/707,451 and Rubinstein et al., “Complex-Valued Eulerian Motion Modulation,” U.S. application Ser. No. 13/607,173, which are all hereinafter incorporated by reference in their entirety. While these methods make small motions visible, the goal of the embodiments of the present invention is to extract quantitative information about heartbeats.

Using Newton's Third Law of Motion to measure cardiac activity dates back to at least the 1930s in using ballistocardiogram (BCG) (Starr, I. et al., “Studies on the estimation of cardiac output in man, and of abnormalities in cardiac function, from the hearts recoil and the bloods impacts; the ballistocardiogram,” The American Journal of Physiology, 127: 1-28 (1939). (Hereinafter, “Starr”)). A subject can be placed on a low-friction platform, and the displacement of the platform due to cardiac activity can be measured. The BCG is not widely used anymore in clinical settings. Other clinical methods using a pneumatic chair and strain sensing foot scale have also been successful under laboratory conditions (Kim, K. et al., “A new method for unconstrained pulse arrival time measurement on a chair,” J. Biomed. Eng. Res., 27: 83-88 (2006). (Hereinafter, “Kim”); Inan, O. et al., “Robust ballistocardiogtam acquisition for home monitoring,” Physiological Measurement, 30: 169-185 (2009). (Hereinafter, “Inan”)). Ballistocardiographic head movement of the sort studied here has generally gained less attention. Such movement has been reported during studies of vestibular activity and as an unwanted artifact during MRI studies (Bonmassar, G. et al., “Motion and ballistocardiogram artifact removal for interleaved recording of EEG and EPS during MRI,” Neurolmage, 16: 1127-1141 (2001). (Hereinafter, “Bonmassar”)). Recently, He, et al. proposed exploiting head motion measured by accelerometers for heart rate monitoring as a proxy for traditional BCG.

An embodiment of the present invention extracts a pulse rate and series of beat sequences from video recordings of head motion. Then, the method evaluates the extracted heart rate and beat location measurements on subjects against an electrocardiogram. Results show that the present method extracts accurate heart rates and can capture clinically relevant variability in cardiac activity.

FIG. 2 is a flow diagram 200 illustrating a high level overview of an embodiment of the present invention. The recorded subject is assembled to be stationary and sitting upright for the duration of the video. First, a region is selected within the head by locating the head region and feature points are tracked for all frames (t=1 . . . T) of the video (202). The method finds a region of interest containing the head and track feature points within the region (e.g., 202). For videos where the front of the face is visible, the method can use the Viola Jones face detector (Viola, P. and M. Jones, “Rapid object detection using a boosted cascade of simple features,” Proceedings of IEEE Conference on Computer Vision and Pattern Recognition, pages 511-518 (2001). (Hereinafter, “Viola”)) from OpenCV 2.4 (Bradski, G., “The Open CV Library,” Dr. Dobb's Journal of Software Tools (2000). (Hereinafter, “Bradski”)) to first find a rectangle containing the face. The method can use the middle 50% of the rectangle widthwise and 90% heightwise from top in order to ensure the entire rectangle is within the facial region. The method also removes the eyes from the region so that blinking artifacts do not affect the results. To do this, the method can remove subrectangle (e.g., a subrectangle spanning 20% to 55% heightwise). For videos where the face is not visible, the method marks the region manually with user input.

The movement of the head throughout the video is measured by selecting and tracking feature points within the region. The method can the OpenCV Lucas Kanade tracker between frame 1 and each frame t=2 . . . T to obtain the location time-series

x_(n)(t) y_(n)(t)

for each point n. Since a modem ECG device operates around 250 Hz to capture heart rate variability and videos are typically shot at 30 Hz, (or frames per second) a cubic spline interpolation can be applied to increase the sampling rate of each y_(n)(t) to 250 Hz.

Many of the feature points can be unstable and have erratic trajectories. To retain the most stable features, the method determines the maximum distance (rounded to the nearest pixel) traveled by each point between consecutive frames, and discards points with a distance exceeding the mode of the distribution.

Both the vertical and horizontal components are extracted from each feature point trajectory (e.g., x₁(t), y₁(t), x_(N)(t), y_(N)(t)) where N is the total number of feature points) (204). In other embodiments, however, the method may extract a vertical component, a horizontal component, or a different component.

Each trajectory is then temporally filtered to remove extraneous frequencies that may be outside the range of possible pulse rates (206). Not all frequencies of the trajectories are required or useful for pulse detection. A normal adult's resting pulse rate falls within the range of [0.75, 2] Hz, or [45, 120] beats/min. Frequencies lower than 0.75 Hz can negatively affect the method's performance, because low-frequency movements like respiration and changes in posture have high amplitude and dominate the trajectories of the feature points. However, harmonics and other frequencies higher than 2 Hz provide useful precision needed for peak detection. Taking these elements into consideration, each y_(n)(t) is filtered to a passband of [0.75, 5] Hz. The method employs a 5^(th) order butterworth filter for its maximally flat passband.

Principal Component Analysis (PCA) decomposes the trajectories into a set of independent source signals (s1(t), s2(t) . . . s5(t)) that describe the main elements of the head motion (208). The underlying source signal of interest is the movement of the head caused by the cardiovascular pulse. The feature point trajectories are a mixture of this movement as well as other motions caused by sources like respiration, vestibular activity and changes in facial expression. The method decomposes the mixed motion into subsignals to isolate pulse by considering the multidimensional position of the head at each frame as a separate data point and using PCA to find a set of main dimensions along which the position varies. The method then selects a dimension on which to project the position time-series to obtain the pulse signal.

Formally, given N feature points, the methods represent the N-dimensional position of the head at frame t as m_(t)=[y₁(t), y₂(t), . . . y_(n)(t)]. The mean and the covariance matrix of the positions are:

$\begin{matrix} {m = {\frac{1}{\lbrack T\rbrack}{\sum\limits_{t = 1}^{T}\; m_{t}}}} & (1) \\ {\sum\limits_{m}\;{= {\frac{1}{T}A{\sum\limits_{t = 1}^{T}\;{\left( {m_{t} - m} \right)\left( {m_{t} - m} \right)^{T}}}}}} & (2) \end{matrix}$

PCA finds the principal axes of variation of the position as the eigenvectors of the covariance matrix: Σ_(m)ϕ_(m)=ϕ_(m)Λ_(m)  (3) where Λ_(m) denotes a diagonal matrix of the eigenvalues λ₁, λ₂, . . . , λ_(N) corresponding to the eigenvectors in the columns of Φ_(m), ϕ₁, ϕ₂, . . . , ϕ_(N).

The method selects the component having clearest main frequency, which identifies the average pulse rate (210). Then, peak detection in the time-domain identifies the beats of the selected signal to calculate beat duration (212).

FIG. 3 displays the two eigenvectors (ϕ₁ and ϕ₂) for two of the subjects (302, 306, and 304, 308, respectively). Each eigenvector, shown as white arrows over the photos of the subjects 302, 304, 306 and 308, represents the N-dimensional direction and magnitude of movement for the feature points. The eigenvectors differ for each subject. The method obtains the one-dimensional (1D) position signal s_(i)(t) by projecting the position time-series onto Ø_(i):

$\begin{matrix} {{s_{i}(t)} = {\begin{pmatrix} m_{1} \\ m_{2} \\ \vdots \\ m_{t} \end{pmatrix} \cdot \phi_{i}}} & (4) \end{matrix}$

During some periods in the video, the head moves abnormally (e.g., swallowing, adjustments in posture). Such movement adds variance to the position vectors, thereby affecting the PCA decomposition. One way to deal with this is discarding a percentage a of the m_(t) with the largest L2-norms before performing PCA. However, all of the m_(t) must still be used in the projection step (Eq. 4) to produce a complete signal. The method sets at 25%.

In an embodiment, an alternative to PCA is independent component analysis (ICA).

The method needs to determine which eigenvector to use for pulse signal extraction. The eigenvectors are ordered such that ϕ₁ represents the most variance in the data, Ø₂ represents the second most, and so on. Although ϕ₁ represents most of the variance, s₁ may not be the clearest pulse signal (e.g., most periodic) for analysis. The method instead chooses the s_(i) that is most periodic. The method quantifies a signal's periodicity as the percentage of total spectral power accounted for by the frequency with maximal power and its first harmonic.

In an embodiment, it is not necessary to consider any signals beyond the first five, i.e., s₁, . . . , s₅ for any of the subjects. The method labels the maximal frequency of the chosen signal f_(pulse) and approximates the pulse rate as 60*f_(pulse) beats per minute.

Average pulse rate alone is not sufficient to fully evaluate the cardiovascular system. Clinicians often assess beat-to-beat variations to form a complete picture. To allow for such analysis, the method performs peak detection on the selected PCA component signal. The peaks are close to

$\frac{1}{f_{pulse}}$ seconds apart with some pulse variability due to the natural variability of heartbeats, variations of the head motion, and noise. The method labels each sample in the signal as a peak if it is the largest value in a window centered at the sample. The method sets the length of the window (in samples) to be round

$\frac{f_{sample}}{f_{pulse}},$ where f_(sample)=250 Hz.

The method can be implemented in MATLAB, for example. Videos can be shot with a Panasonic Lumix GF2 camera in natural, unisolated environments with varying lighting. The videos can have a frame rate of 30 frames per second, 1280×720 pixel resolution and a duration of 70-90 seconds in one embodiment. Subjects can be connected to a wearable ECG monitor (Delano, M., “A long term wearable electrocardiogram (ECG) measurement system,” Master's Thesis, Massachusetts Institute of Technology (2012) (hereinafter, “Delano”)) to compare the results of the method to the ECG, for testing purposes. This device has a sampling rate of 250 Hz and three electrodes which are placed on the forearms, in one embodiment.

The method extracted pulse signals from 18 subjects with a frontal view of the face (as in FIG. 3). The subjects varied in gender (7 female, 11 male) and skin color. They ranged from 23-32 years of age and were all seemingly healthy. The method calculated the average pulse rate using the frequency of maximal power for the selected PCA component. Similarly, the method computed the true pulse rate by finding the main frequency of the ECG spectrum. Table I presents the results. The average rates are nearly identical to the true rates for all subjects, with a mean error of 1.5%. The number of peaks were also close to ground truth values, with a mean error of 3.4%.

The method evaluates the ability of the signal to capture subtle heart rate variability. Clinically meaningful HRV measures typically use 10-24 hours of ECG data, therefore testing of the method did not attempt to compute any of these for the 90 second videos. Instead, the distributions of time between successive peaks for each signal are compared for testing. Incorrect or missed peaks can introduce spurious intervals too large or small to be caused by the natural variations of the heart. Therefore, only intervals with length within 25% of the average detected pulse period are considered.

The Kolmogorov-Smirnov (KS) test measures the similarity of the distributions, with the null hypothesis being that the observations are from the same distribution. Table 2 presents the results. At a 5% significance level, 16 of the 18 pairs of distributions were found to be similar.

FIG. 4 is a diagram 400 illustrating histograms of four of the 16 distributions binned at every 0.05 seconds. The method captured a wide range of beat-length distributions shapes, from the flat distribution of subject 4 (402, 404) to the peakier distribution of subject 10 (414, 416). FIG. 4 shows each subject having an ECG histogram (402, 406, 410, 414) being similar to a respective motion based histogram (404, 408, 412, 416) of an embodiment of the present invention

TABLE 1 Average pulse rate and # peaks detected from ECG and by the method. Avg. Pulse (beats per minute) Number of beats Sub. ECG Motion (% error) ECG Motion (% error) 1 66.0 66.0 (0)   99 98 (1.0) 2 54.7 55.3 (1.1) 82 84 (2.4) 3 81.3 82.6 (1.6) 122 116 (4.9)  4 44.7 46.0 (2.9) 67 70 (4.5) 5 95.3 96.0 (0.7) 143 142 (0.7)  6 78.9 78.0 (1.1) 92  78 (15.2) 7 73.3 71.3 (2.7) 110 100 (9.1)  8 59.3 58.6 (1.2) 89 88 (1.1) 9 56.7 58.6 (3.4) 85 84 (1.2) 10 78.7 79.3 (0.8) 118 117 (0.8)  11 84.7 86.6 (2.2) 127 121 (4.7)  12 63.3 62.6 (1.1) 95 95 (0)   13 59.3 60.0 (1.2) 89 89 (0)   14 60.0 61.3 (2.2) 90 89 (1.1) 15 80.0 81.3 (1.6) 120 114 (5.0)  16 74.7 74.6 (0.1) 112 110 (1.8)  17 50.0 49.3 (1.4) 75 76 (1.3) 18 77.1 78.8 (2.2) 90 85 (5.6)

Table 2 presents results when comparing the interpeak distributions of the ECG and the method including the means (μ) and standard deviations (σ) of each distribution, the number of outliers removed from the distribution, and the p-value of distribution similarity. 16 of the 18 pairs of distributions were not found to be significantly different.

TABLE 2 ECG Motion KS-Test Sub. μ(σ) μ(σ) p-value 1 0.91(0.06) 0.90(0.06) 0.89 2 1.08(0.08) 1.06(0.11) 0.52 3 0.73(0.04) 0.73(0.08) 0.05 4 1.34(0.19) 1.28(0.18) 0.14 5 0.62(0.03) 0.63(0.07) <0.01 6 0.76(0.04) 0.76(0.04) 0.64 7 0.81(0.05) 0.81(0.06) 0.85 8 1.01(0.04) 1.02(0.09) 0.16 9 1.04(0.07) 1.04(0.11) 0.27 10 0.75(0.04) 0.75(0.04) 0.75 11 0.70(0.06) 0.70(0.08) 0.30 12 0.94(0.08) 0.94(0.09) 0.85 13 0.99(0.04) 0.98(0.12) <0.01 14 0.99(0.11) 0.98(0.12) 0.47 15 0.74(0.05) 0.75(0.06) 0.95 16 0.80(0.05) 0.80(0.06) 0.60 17 1.18(0.08) 1.18(0.11) 0.70 18 0.76(0.05) 0.76(0.06) 0.24

Pulse motion constitutes only a part of total involuntary head movement. The magnitude of the different movements within [0.75, 5] Hz is quantified by calculating root mean square (RMS) amplitudes of the feature point trajectories. For each subject, the mean RMS amplitude of the trajectories is calculated before and after filtering to a passband within 5% of the pulse frequency. The mean RMS amplitude of the trajectories without filtering was 0.27 (std. dev of 0.07) pixels across the subjects. The mean RMS amplitude after filtering to the pulse frequency was 0.11 (0.05) pixels. Thus the pulse motion had roughly 40% the RMS amplitude of other head motions within the [0.75, 5] Hz frequency range. The robustness of the method can be compared to a color-based pulse detection system (Poh, M. et al.) in the presence of noise. The color method spatially averages the R, G, and B channels in the facial area and uses independent component analysis (ICA) to decompose the signals into three independent source signals. The source with the largest peak in the power spectrum is then chosen as the pulse signal.

Varying levels of zero-mean Gaussian noise can be added to the videos and swept the standard deviation from 5 to 500 pixels. For each subject, σ_(motion), the maximum noise standard deviation, was calculated before the method first produced an average pulse rate outside 5% of the true rate. σ_(color) was calculated in a similar manner for the color method.

FIG. 5 is a diagram 500 illustrating the results of the subjects, plotting σ_(color) based on σ_(motion). Each point represents one of the results of the subjects σ_(color) value based on σ_(motion). The method outperformed color for seven of the seventeen subjects, and performed worse for nine subjects. Note that color failed to produce a correct pulse rate for subject seven before adding any noise.

There is a large variance in σ_(motion) and σ_(color) across the subjects, suggesting that there are subject-specific factors that affect performance. To understand why, σ_(motion) is compared to β, the ratio of the total energy of the feature points within 5% of f_(pulse) to the maximal energy at any other frequency.

FIG. 6 is a graph 600 illustrating a plot of σ_(motion) against β for all subjects on a logarithmic scale. The subjects with the 10 highest σ_(motion) values also have 10 of the top 11 β values. This indicates that the method performs best for subjects with a large ballistocardiac motion relative to any other periodic head movement.

No similar relationship between σ_(color) and the frequency content of the Red, Green, and Blue (RGB) channels was found, likely because of the layer of complexity introduced by the ICA algorithm. However, when simplifying the method to extracting a signal from the G channel alone, noise performance is shown to be strongly related to the ratio of power at the pulse frequency to the next largest power in the spectrum. No relationship was found between motion or color performance and skin tone.

FIG. 7 are frames 702, 704, and 706 from videos employed by an embodiment of the present invention. One advantage of motion-based detection over color is that a direct view of the skin is not needed. Videos of the backs of the heads (702, 704) of eleven subjects and a video of one subject wearing a mask (706), as shown in FIG. 7, were able to produce average heart rates close to the true values for all videos.

FIG. 8 is a diagram showing a frame 802 from a 30-second video of a newborn recorded in situ at a hospital nursery and corresponding graphs 804 and 806. A video of the baby's actual pulse rate is shown on a hospital-grade monitor 808 measuring the perfusion of blood to its skin. The method extracts a clean pulse signal 804 and a corresponding frequency response that matches the pulse rate reported by the monitor.

Embodiments of the present invention consistently obtain accurate pulse rate measurements from head motion. The results for beat detection are equally encouraging. Most beat interval distributions looked qualitatively similar to the ECG distributions, indicating that the method captures a real physiological variability. For sixteen of the eighteen subjects, there is no statistically significant difference between the ECG and the motion beat intervals. This is a stronger test than is required in most clinical contexts. Typically heart rate variability (HRV) is used to dichotomize patients into high and low risk groups, so the precise shape of the distribution is not relevant. The relevant test is whether the distribution of motion-generated intervals yields the same set of high risk individuals as ECG generated intervals. Since all subjects were healthy volunteers, high risk individuals were not tested.

Several factors affected the results. First, the camera has a sampling rate of 30 Hz. ECG used for HRV analysis normally has a sampling rate of at least 128 Hz. Cubic interpolation of the signal only partially addresses this discrepancy. Second, extra variability might be introduced during the pulse transit time from the abdominal aorta to the head. In particular, arterial compliance and head mechanics could affect the results. Third, the variable and suboptimal lighting conditions can affect the feature tracking Finally, the videos were only a maximum of 90 seconds long. Normally, HRV measures are computed over many hours to obtain reliable estimates.

An important future direction is to develop approaches for moving subjects. This is complicated because, as the results show, even other involuntary head movements are quite large in relation to pulse motion. Clearly with larger motions such as talking, more sophisticated filtering and decomposition methods are needed to isolate pulse.

This method considers the frequency and variability of the pulse signal. However, head movement can offer other information about the cardiac cycle. If head displacement is proportional to the force of blood being pumped by the heart, it may serve as a useful metric to estimate blood stroke volume and cardiac output. Additionally, the direction of the movement could reveal asymmetries in blood flow into or out of the head. This might be useful for the diagnosis of a stenosis, or blockage, of the carotid arteries.

Another future direction is to better assess the strengths and weaknesses of the color and motion pulse estimation methods. The results suggest that neither method is strictly more robust than the other in the presence of noise. However, further work needs to be done with varying lighting, skin tones, and distance from the camera to form a complete picture. In addition, sensitivity of the method to voluntary motions like talking or typing. For many applications, this is a critical factor. A motion-based approach is certainly better when the face is not visible. Based on these ideas, a combination of the color and motion methods are more useful and robust than using either one independently.

The present invention offers a non-invasive, non-contact means of cardiac monitoring. The method takes video as input and uses feature tracking to extract heart rate and beat measurements from the subtle head motion caused by the Newtonian reaction to the pumping of blood at each heartbeat. A combination of frequency filtering and PCA allows identification of the component of motion corresponding to the pulse and then extraction of peaks of the trajectory to identify individual beats. When evaluated on eighteen subjects, the method produced virtually identical heart rates to an ECG and captured some characteristics of inter-beat variability.

FIG. 9 is a block diagram 900 illustrating an example embodiment of the present invention. A subject 902 is videotaped by a camera 904 for a period of time, recording a video 906 having the subject's 902 head in the frame. A selection module 908 receives the video 906 recorded by the camera 904 directly or via in intermediary. The selection module 908 selects a region of the video 910 having the subject's head and forwards it to a tracking module 912. The selection of the region of the video 910 is described above. The tracking module 912 then tracks features of the video 914 which are forwarded to the analysis module 916. The analysis module calculates an oscillation or heart rate 918 based on the features 914, which it outputs to a display 920.

Embodiments or aspects of the present invention may be implemented in the form of hardware, software, or firmware. If implemented in software, the software may be any form of software capable of performing operations consistent with the example embodiments disclosed herein. The software may be stored in any non-transitory computer readable medium, such as RAM, ROM, magnetic disk, or optical disk. When loaded and executed by processor(s), the processor(s) are configured to perform operations consistent with the example embodiments disclosed herein. The processor(s) may be any form of processor(s) capable of being configured to execute operations as disclosed herein.

The teachings of all patents, published applications and references cited herein are incorporated by reference in their entirety.

While this invention has been particularly shown and described with references to example embodiments thereof, it will be understood by those skilled in the art that various changes in form and details may be made therein without departing from the scope of the invention encompassed by the appended claims. 

What is claimed is:
 1. A method comprising: tracking trajectories of respective feature points within a selected region of a digital video by: determining a maximum distance of movement of each feature point between consecutive frames of the digital video; and discarding feature points corresponding to points having a maximum distance between consecutive frames of the digital video exceeding a mode of distribution; extracting a plurality of dimensional components from each tracked trajectory of the respective feature points within the selected region of the digital video; and determining an oscillation rate of a subject represented by the feature points of the digital video by analyzing the dimensional components of the tracked trajectories of the respective feature points within the selected region of the digital video.
 2. The method of claim 1, wherein the oscillation rate is a heart rate or pulse rate of the subject.
 3. The method of claim 1, wherein the plurality of extracted components includes a vertical component, horizontal component, or vertical and horizontal components.
 4. The method of claim 1, wherein selecting the region of the digital video includes finding a region of interest including a head of the subject.
 5. The method of claim 1, wherein analyzing the motion of the feature points includes: filtering the extracted feature points temporally; generating a plurality of components using principal component analysis on the temporally filtered extracted feature points; selecting one of the components of the plurality of components generated by the principal component analysis; detecting peaks of the selected component; and determining the oscillation rate based on the detected peaks of the selected component.
 6. The method of claim 5, wherein filtering the extracted feature points temporally removes frequencies that are outside a predetermined range of pulse rates.
 7. The method of claim 5, wherein performing principal component analysis includes decomposing the filtered feature points into a set of independent source signals.
 8. The method of claim 5, wherein selecting the component includes selecting the component having a clearest main frequency.
 9. A system comprising: a processor; and a memory with computer code instructions stored therein, the memory operatively coupled to said processor such that the computer code instructions configure the processor to implement: a tracking module configured to track trajectories of respective feature points within a selected region of a digital video by: determining a maximum distance of movement of each feature point between consecutive frames of the digital video; and discarding feature points corresponding to points having a maximum distance between consecutive frames of the digital video exceeding a mode of distribution; an extraction module configured to extract a plurality of dimensional components from each tracked trajectory of the respective feature points within the selected region of the digital video; and an analysis module configured to determine an oscillation rate of a subject represented by the feature points of the digital video by analyzing the dimensional components of the tracked trajectories of the respective feature points within the selected region of the video.
 10. The system of claim 9, wherein the oscillation rate is a heart rate or pulse rate of the subject.
 11. The system of claim 9, wherein the extracted component is a vertical component, horizontal component, or vertical and horizontal component.
 12. The system of claim 9, wherein the selected region is selected by a selection module implemented by the processor, and the selection module is further configured to find a region of interest including a head of the subject.
 13. The system of claim 9, wherein the analysis module is further configured to filter the extracted feature points temporally; generate a plurality of components using principal component analysis on the temporally filtered extracted feature points; select one of the components of the plurality of components generated by the principal component analysis; detect peaks of the selected component; and determine the oscillation rate based on the detected peaks of the selected component.
 14. The system of claim 13, wherein the analysis module is further configured to filter the extracted feature points temporally by removing frequencies that are outside a predetermined range of pulse rates.
 15. The system of claim 13, wherein the analysis module is further configured to perform principal component analysis by decomposing the filtered feature points into a set of independent source signals.
 16. The system of claim 13, wherein the analysis module is further configured to select the component by selecting the component having a clearest main frequency.
 17. A non-transitory computer readable medium having instructions stored thereon, the instructions configured to, when loaded and executed by a processor, cause the processor to: track trajectories of respective feature points within a selected region of a digital video by: determining a maximum distance of movement of each feature point between consecutive frames of the digital video; and discarding feature points corresponding to points having a maximum distance between consecutive frames of the digital video exceeding a mode of distribution; extract a plurality of dimensional components from each tracked trajectory of the respective feature points within the selected region of the digital video; and determine an oscillation rate of a subject represented by the feature points of the digital video by analyzing the dimensional components of the tracked trajectories of the respective feature points within the selected region of the digital video. 